DALS006-统计推断(Inference)05-置换检验(Permutation Test)

发表于 | 分类于
| 字数统计: 2,577 | 阅读时长 ≈ 10

前言

这篇笔记是《Data Analysis for the Life Sciences》的第2章：统计推断(Inference)的第5部分。这一部分的主要内容涉及置换检验。

置换检验英文是Permutation test，这是Fisher于20世纪30年代提出的一种基于大量计算（computationally intensive），利用样本数据的全（或随机）排列，进行统计推断的方法，因其对总体分布自由，应用较为广泛，特别适用于总体分布未知的小样本资料，以及某些难以用常规方法分析资料的假设检验问题。在具体使用上它和Bootstrap Methods类似，通过对样本进行顺序上的置换，重新计算统计检验量，构造经验分布，然后在此基础上求出P-value，进行推断，生信中常用的GSEA分析本质上就是置换检验。

置换检验的根本原理就是重新采样增加小样本的整体样本量，然后看其的概率分布来预测假设是否成立。

置换检验(Permutation Test)

现在考虑一种情况，在这种情况下，标准的数理统计近拟(standard mathematical statistical approximation)都不适用。但是，我们已经计算出了一个总结性的统计，例如均值的差值，但是，我们的这些数据不符合一些即定分布，也没有办法得到总体，也就是说无法进行一些数据模拟，例如Monte Carlo模拟，面对这种情况，我们可以采用置换检验(Permutation test)。

置换检验过程

在前面部分中，我们介绍了参数模拟，这种方法有助于我们计算出观察到的差值是否有统计学的意义。而置换检验则会使用另外一种计算方法，它也是利用了两组的差值。在这种方法中，我们会随机置换干预组与对照组的分组标签。当我们置换了干预组与对照组的分级标签后，得到的新数据我们会假设它们服从零分布，现在我们通过置换这两组数据来计算一下零分布（这里假设置换1000次）。

第一步，是先计算出两组（这里是treatment组与control组） 的均值差值，如下所示：

1
2
3
4
5
6
7
library(dplyr)
dir <- system.file(package = "dagdata")
filename <- file.path(dir,"extdata/femaleMiceWeights.csv")
dat <- read.csv(filename)
control <- filter(dat, Diet == "chow") %>% dplyr::select(Bodyweight) %>% unlist
treatment <- filter(dat, Diet == "hf") %>% dplyr::select(Bodyweight) %>% unlist
obsdiff <- mean(treatment) - mean(control)

第二步，对样本重新进行分组，也就是将上述的两组的所有数据（每组都是12个样本，一共是24个数据）重新进行分组，每组12个数值，每组重新分组时，都计算出均数差，24个数据随机等分为2组，一共是2704156个均数差（使用choose(24,12)就能计算出来，不过我们这里只不进行全部置换，我们只置换1000次，现在我们通过置换这两组数据来计算一下零分布（置换1000次），如下所示：

1
2
3
4
5
6
7
8
9
N <- 12
avgdiff <- replicate(1000, {
  all <- sample(c(control, treatment))
  newcontrols <- all[1:N]
  newtreatments <- all[(N+1):(2*N)]
  return(mean(newtreatments)-mean(newcontrols))
})
hist(avgdiff)
abline(v=obsdiff,col="red",lwd=2)

第三步：计算出上面那1000个均数差中大于obsdiff的个数。

现在我们考虑一下，有多少个零均值(null means，结合上下文，这里的零均值应该是在零假设成立的前面下的均值）大于观察值呢？这个比值就是零假设成立下的p值。我们向分子与分母都加1，用于校正一下p值（书中列出的参考文献，这里不解释它的数学原理），如下所示：

1
2
# the proportion of permutations with larger difference
(sum(abs(avgdiff) > abs(obsdiff))+1)/(length(avgdiff)+1)

计算结果如下所示：

1
2
3
> # the proportion of permutations with larger difference
> (sum(abs(avgdiff) > abs(obsdiff))+1)/(length(avgdiff)+1)
[1] 0.04795205

第四步：就是根据第三步中的计算结果（上面的计算结果是0.04795205）来下结论。

如果p>0.05，说明在两样本来自同一总体的假设下，当前样本的出现是很平常的，不能拒绝H0假设，认为两样本的差异无统计学意义；如果p<0.05，则认为两样本的差异有统计学意义，不过考虑到这里的p值基本上就是在0.05附近，因此下结论也应该慎重。

现在我们使用小样本数据来重复一下这个实验，我们通过抽样的方式来创建一个小样本数据，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
N <- 5
control <- sample(control, N)
treatment <- sample(treatment, N)
obsdiff <- mean(treatment) - mean(control)
avgdiff <- replicate(1000,{
  all <- sample(c(control, treatment))
  newcontrols <- all[1:N]
  newtreatments <- all[(N+1):(2*N)]
  return(mean(newtreatments) - mean(newcontrols))
})
hist(avgdiff)
abline(v=obsdiff, col="red",lwd=2)

上面的直方图显示的是使用置换检验计算出来的小样本数据的均值之差，红线表示的是观察值的差异。

现在我们发现，使用这种方法得到的观察差值并不显著，这里需要记住，没有理论能保证，置换的结果接近于实际的零分布。例如，如果不同的总体之间确实存在着差异，一些置换检验的结果可能就计算不出差异，也有一些能计算出差异。这说明，在置换检验中，零分布比实际的零分布有着一个更扁平的尾部。这也就说明了，为什么置换检验会得到一个比较保守的p值。出于这种原因，当我们的样本数目比较少时通常不做置换检验。

另外我们还要注意，置换检验也有是假设的：假设样本是独立的，是可以进行交换的(“exchangeable”)。如果你的数据中含有隐藏的结构，置换检验生成的估计的零分布分则会低估该分布的尾部面积，这有可能破坏原始数据中的现在结构（注：这里我的理解就是，假设A样本和B样本数据没有差异，那么它们无论怎么置换后，得到的数据应该是差不多的，也就是说得到的p值会比较大）。

练习

书中还列出一个练习题，现在我们来做一下。

第一步，下载原始数据，如下所示：

1
2
3
4
5
6
7
8
url <- "https://raw.githubusercontent.com/genomicsclass/dagdata/master/inst/extdata/babies.txt"
filename <- basename(url)
library(downloader)
library(dplyr)
download(url, destfile=filename)
babies <- read.table("babies.txt", header=TRUE)
bwt.nonsmoke <- filter(babies, smoke==0) %>% dplyr::select(bwt) %>% unlist
bwt.smoke <- filter(babies, smoke==1) %>% dplyr::select(bwt) %>% unlist

第二步，现在我们生成一个数据集，先来观察一下它们的差值，如下所示：

1
2
3
4
5
6
N=10
set.seed(1)
nonsmokers <- sample(bwt.nonsmoke , N)
smokers <- sample(bwt.smoke , N)
obs <- mean(smokers) - mean(nonsmokers)
obs

结果如下所示：

1
2
> obs
[1] -16.9

现在我们的问题就是：在我们构建的这个数据里（这个数据含有1个nonsmokers样本，1个smokers样本），这个观察到的差值-16.9是否有统计学上的意义？

这里我们使用前面提到的置换检测来计算一下上面构建好的p值。我们会重新杂揉(原文是reshuffle，没找到相应的中文翻译，这里我译为“杂揉”，也就是把这两组数据混合起来)，然后重新计算它们的均值。现在我们称使用下面的代码来计算一下，经过置换后的，其中的一个均值的差值，如下所示：

1
2
3
4
5
dat <- c(smokers,nonsmokers)
shuffle <- sample(dat)
smokersstar <- shuffle[1:N]
nonsmokersstar <- shuffle[(N+1):(2*N)]
mean(smokersstar)-mean(nonsmokersstar)

结果如下所示：

1
2
> mean(smokersstar)-mean(nonsmokersstar)
[1] -8.5

上面的这个值就是我们即将构建的零分布中的一个观测值。

现在我们开始计算p值，将随机数种子设为1，即set.seed(1)，用于置换1000次，生成一个零分布。现在我们从这个零分布中观察到的前面的那个差值的p值是多少呢？

计算过程的完整代码如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
N=10
set.seed(1)
nonsmokers <- sample(bwt.nonsmoke , N)
smokers <- sample(bwt.smoke , N)
obs <- mean(smokers) - mean(nonsmokers)
avgdiff <- replicate(1000,{
  shuffle <- sample(dat)
  sample_smokers <- shuffle[1:N]
  sample_nonsmokers <- shuffle[(N+1):(2*N)]
  return(mean(sample_smokers) - mean(sample_nonsmokers))
})
hist(avgdiff)
abline(v=obs, col="red",lwd=2)
(sum(abs(avgdiff) > abs(obs))+1)/(length(avgdiff)+1)

结果如下所示：

1
2
> (sum(abs(avgdiff) > abs(obs))+1)/(length(avgdiff)+1)
[1] 0.05394605

从上面我们可以看到，p值是大于0.05的。我们来看一下直方图的结果：

1
2
hist(avgdiff)
abline(v=obs, col="red",lwd=2)

现在我们分别使用t检验和Wilcoxon检验分别计算如下，如下所示：

1
2
t.test(smokers,nonsmokers)
wilcox.test(smokers,nonsmokersstar)

结果如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
> t.test(smokers,nonsmokers)
	Welch Two Sample t-test
data:  smokers and nonsmokers
t = -2.067, df = 13.159, p-value = 0.059
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -34.5419234   0.7419234
sample estimates:
mean of x mean of y 
    109.9     126.8 
> wilcox.test(smokers,nonsmokersstar)
	Wilcoxon rank sum test with continuity correction
data:  smokers and nonsmokersstar
W = 34, p-value = 0.2399
alternative hypothesis: true location shift is not equal to 0

从上面的结果来看，置换检验算出的p值为0.053，t检验计算出的p值为0.059，这两个计算结果还是比较接近的。最后我们使用了Wilcoxon检验，p值是大于0.05的，毕竟Wilcoxon检验的功效较小。

练习中的第2题提到，我们不使用均值之差，而是使用中位数之差，即obsdiff <- median(smokers) - median(nonsmokers)，看一下p值，计算过程如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
N=10
set.seed(1)
nonsmokers <- sample(bwt.nonsmoke , N)
smokers <- sample(bwt.smoke , N)
(obs <- median(smokers) - median(nonsmokers))
dat <- c(nonsmokers,smokers)
  
avgdiff <- replicate(1000, {
  shuffle <- sample(dat)
  smokersstar <- shuffle[1:N]
  nonsmokersstar <- shuffle[(N+1):(2*N)]
  obs = median(smokersstar)-median(nonsmokersstar)
  return(obs)
})
hist(avgdiff)
abline(v=obs, col="red",lwd=2)
(sum(abs(avgdiff) > abs(obs))+1)/(length(avgdiff)+1)

结果如下所示：

1
2
> (sum(abs(avgdiff) > abs(obs))+1)/(length(avgdiff)+1)
[1] 0.02197802

从结果来看，通过来看中位数差值，p值是小于0.05的。来看一下直方图的结果：

出现这种结果我的理解就是因为，随机了取样的数目太少，是10个，如果变成20个，有可能p值就大于0.05了（当把N变为20，那么p值就在于0.05，并且中位数的直方图也更像正态分布的钟形曲线）。

参考资料

  1. Permutation Test 置换检验(转）